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ABSTRACT 

The problem of maximum range atmospheric reentry for an orbiting 
lifting glider was treated by Bryson and Denhan “Gy the method of 
"steepest ascent". The same problem is undertaken here by a method 
of differential corrections developed by Faulkner. This method makes 
use of a Newton-Raphson type iteration based on paths which satisfy 
the Euler-Lagrange equations. A comparison of results is made, show- 
ing large differences in control variable history, and longer range 
for the path obtained by differential corrections. The problem was 
characterized by a sharp "ridge" in the domain of the starting values 


of the adjoint variables and the effect of this on the convergence of 


both methods is discussed. Finally, the difficulty of choosing initial 


approximations for the starting values of the adjoint variables is 


discussed, and a method is presented for obtaining these from a nominal 


path as the first step in the computer routine. 


iL 


A. E. Bryson and W. F. Denham, A steepest-ascent method for solve 


ing optimum programming problems (Ratheon report BR-1303), Raytheon 
Company, Bedford, Mass., 1961. 
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INTRODUCTION 


A problem of optimum control may be described very roughly as 
the determination of the decision=making or control function which 
Will result in the largest possible value of the final "payoff" from 
some continuing process, 

The process might be illustrated as shown in figure 1.1, where 
P is the "payoff", something whose value depends on the terminal point 


or the path history, and Q is a set of given constraints which must be 


Satisfied. 

a 
ra 
|— 

4 Given initial point eee 
aon, — \ ce 
| x f End point a 
ye t = ib 
P = maximum 
—_ Q=0 


Fig. 1.1 An Optimum Path 


The methods of solution commonly used now were foreseen at least 
forty years ago, but the volume of computation required limited their 
application to only the simplest of problems. In the past ten years 
there has been an explosion of interest due to the availability of 
highespeed digital computers on the one hand, and on the other hand 
to the urgent need for solutions in such applications as space 


trajectories. The "state of the art" is such that no method is com- 








pletely general, and a great diversity of approaches is employed, 
reflecting the absence of a firm framework of mathematical theory 
concerning nonlinear differential equations. 

in this thesis a comparison 1s made between two methods for cal- 
culating optimum angle of attack programs for a lifting unpowered 
vehicle to attain maximum range, starting from circular orbital 
velocity at a height of 300 thousand feet above the surface of the 
earth, The problem was chosen because it is typical of a class of 
problems in which there is now great interest, and because two of 
the foremost workers in the field of optimization, Arthur E. Bryson 
and W. F. Denham, had published a clear and well-documented solution 
by the method of "steepest ascent", 

The comparison was Suggested by Stanley Ross, while one of the 
writers was engaged in summer field work under the sponsorship of 
the Lockheed Missile and Space Company. John V. Breakwell and George 
Leitmann generously gave many hours of their time to general dis= 
cussions and made specific suggestions which were very helpful. 
Professor Frank D. Faulkner of the U. S. Naval Postgraduate School 
not only laid out the general form of the problem and overcame the 
difficulties that arose, but also the method of solution used was 
the one developed by him from the basic method of differentials of 


Bliss, 








Chapter I 


GENERAL FORMULATION OF A PROBLEM OF OPTIMUM CONTROL 


1.1 The equations of motion, the control variables and the cone 
straints. 

It will be helpful to define an optimum control problem in 
somewhat more exact terms than those given in the introductory 
remarks, 

Suppose we are given a system whose behavior is described 
by a system of differential equations which we shall call the 
equations of motion. We will suppose that wherever any of these 
differential equations is of higher than first order, we have de- 
fined such additional variables as necessary so as to obtain n 


first order differential equations of the form 


S; = £4(Sqz, Sos ecer Sys Pps Poo coos Duo t) 
T= 1, 25 ooo, N 
where the variables S, are the state variables, the Pj are the 
control variables, t is the independent variable, and the dot 
denotes differentiation with respect to the independent variable. 
It is assumed that the f, are of class 0". 
It will be convenient for what is to follow, if matrix 

notation is adopted at this point. The equations of motion may 
be written as the single matrix equation 


fake) S = F(S, P, t) 
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where 


=a! 
5S =]. |, an nxl matrix of the dependent state 
A variables 
Ss 
n 
Py 
F=], » an nxl matrix of known functions of the 
. control variables, the state variables, 
i and possibly the independent variable 
Py 
Pia as » an mxl matrix of control variables, 
; which we are free to choose as functions 


Pa of t. 


In addition, there are k constraints or end conditions, of the form 
e, (S,t) = 0, which may define values of the state variables at the 


end point t =T, or may be of integral form such as 


aes 
/r, fs, dt = constant. 


In the latter case, we shall define an additional variable s, 4) 


so that 


With the initial and end conditions 





S41 CO) e=a6 


“n+l (T) = the given constant. 


These may be written as 


e 


a 
wr?) E=!. |, a kxl matrix of constraint fimctions. 
ek 


We wish to discover the particular control variable 
program which will drive this system from a given initial point 
to some terminal point which, of the set of all such points which 
Satisfy the terminal constraints, gives a maximum value to some 
function, M(S,t), of the state variables; that is, for the ter- 
minal point t =T, possibly unknown, 


M(S,t)p =M(T) = max. 


1.2 Bounded control variables and inequality constraints. 
The permissible domain of the control variables may likewise 
be limited. Some examples of restrictions on the control for com- 


mon physical systems might be: 


2 Fy SP max, 


a 
{? dt <I max, 
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ee ae Lee 0 et <tj and 
P = 0 for t,<t< ue 
Bounded control variables appear frequently in real physical 
systems. For solving these problems it may be useful to re- 
place the inequalities by mathematically equivalent equalities. 
This is accomplished by the introduction of suitable new vari- 
ables in the manner due to Valentine (ref.1). For instance, 
for a constraint 
OSp<sP max 
one may define the real variable g, so that 
pibieeyaip eae = 0: 
This procedure transforms the inequality to an equality, and we 
adjoin this equation to the set (1.1). The same procedure is 
applied in the case of inequality constraints or end conditions. 

We will generally require that all state variables be speci- 
fied initially, but this is not necessary. One or more of these 
may be unspecified, that is "free", 

The general form of these problems, then, is the two-point 
boundary value problem. With the stipulation that the functional 
M is to be maximized (or minimized) they become, depending on the 
formulation, the classical problem of Bolza or problem of Mayer 
of the calculus of variations. Since to maximize some quantity, 
say u, is equivalent to minimizing the quantity minus u, we will 


henceforth consider that the term "maximum" includes both cases 








unless otherwise specified. 
1.3 Degenerate problems. 

Some important obServations may be made here. First, as 
pointed out by Faulkner (ref.2) and Breakwell (ref.3), the case 
in which all of the functions f, involve any function of a par- 
ticular control variable only linearly are said to be degenerate. 
These degenerate problems typically are of the "bang-bang" type 
in which the control changes discontinuously. The neglect of 
induced drag in aerodynamic problems for instance often results 
in bang-bang control. 

second, the constraints must not involve the control variables, 
or the problem is overspecified. If the control program is to be 
such as to produce an extreme of the functional M, then it may not 
Simultaneously satisfy any prescribed constraint at any isolated 
point of the path, else the solution is generally discontinuous. 
1.4 The adjoint variables. 

We shall introduce a set of Lagrange multipliers, as yet un- 
determined 

U = Ee Uns sees u,,| 

a lx matrix of adjoint variables. The term "adjoint" is used, 
Since they will be chosen to be solutions to the system of differ. 
ential equations which is adjoint to equations for the variations 
of the given system. 


Suppose we write (1.1) as 








and integrate its product with the adjoint vector between the 
limits of the independent variable. For convenience, we will 


call these limits zero and T, and we have, 


Hh @ 
J RUS SPE er). 
© 
The corresponding variational equation is 
eS 
J us - F.6S - F &) dt = 0. 
O Ss p 
Now this is integrated by parts to eliminate from the integrand 


the variations of the state variables. 
(1.3) [us| : =f U+UF_)§S +U F sp | at-Po F| ve at 
O 0 Ss p “k tie ik 


where F, and Bo are the matrices of partial derivatives, 


af, /as, oe af, /2s, 
(1.4) ee ; » an nxn matrix, 


af, /as, . - - of,/as, 


af, /ap, . » » af,/aP_ 


bry 
1] 


Joan MoaiMaeriog, 


af /aP, + + + 3f,/ar, 


and te is a symbol for any point or set of points where S or U 
is discontinuous. 


Now, if U is chosen as a solution to the equation 








(led) U+UF_ = 0, 
then (3) becomes 


es | iE : ilar 
(1.6) [vss] - J (U F.éP) at - © [ur] katy. 


This sequence of operations defines the system which is adjoint 
to the variational equations of (1.1), $s . Fo dS = FOP ; 


and forms the 1 xn matrix, or vector 


U = E Us —_—__ 7 é 


Hereafter, no distinction will be drawn between a vector in 
kespace anda lxk or kx1 matrix, or a corresponding 
column or row of a matrix. 

1.5 Green's formula and the Euler-Lagrange equations. 

Equation (1.6)is the fundamental formula relating the varia- 
tions of the end values of the state variables with the variations 
of the control variables. Note that the formula shows the ter. 
minal value of the adjoint variable to be the sensitivity of the 
terminal value of the corresponding state variable to a change in 
control. This interpretation will be discussed further in con- 
nection with the transversal conditions (sect. l. 9 ). 

Equation (1.6) is often called Green's formula (one form). 


See Coddington and Levinson (ref. 4+) page 86. The last term on 








the right applies only at points of discontinuity of F, and 

Since these are discontinuities in S, they are ®corners", 
Hereafter, it will be assumed that the solution does not have 

a corner unless it is otherwise mentioned. In this case the 

last term drops out. The fundamental lemma of the calculus of 
variations states that, given a control program P(t) such that 
the constraints E are satisfied, then if the coefficient of 6P 

in (1.6) does not vanish for some solution U to the adjoint 

it is possible to construct a variation such that the constraints 
are still satisfied, but the function M exceeds the value on the 
first path. A proof of this lemma is included in refs. (6) and (10). 
Thus for M to have at least a stationary value M*, 

ele) UF, = 0, (The Euler equation) 

Equation (1.5) defines the system of equations which is adjoint 
to the variations of the equations of motion, (1.1). The adjoint 
equations and equation (1.7) will be called the Euler-Lagrange 
equations. It will be shown that these, together with the con- 
straints and the equations resulting from the transversal condi- 


tions comprise a complete set which determines the optimum solution. 


1.6 The maximum principle of Pontryagin. 
The functions fs are often composed of terms which are funce 
tions of the state variables alone and terms which contain also 


the control variables. Suppose we collect terms of the first type 


10 








into ann x 1 matrix we will call V, and those of the second type 
into ann x 1 matrix G. Then (1.1) becomes 

S=V+G. 
Now consider these matrices to be vectors in an n-dimensional 
hyperspace of the state variables, and a particularly fruitful 


geometric interpretation of (1.7) appears. Let us consider equa- 


tion (1.7) 
= alia. 
10 ae ae 
p 
as the condition that 
_"s oes, 
U°*G = max 


p 


that is, let us require that the vector G be everywhere so dir- 
ected that the projection of G upon U is maximized. The control 
program which maintains this relationship (while satisfying the 
constraints) is optimum at least locally. This statement of the 
condition is due to the Russian mathematician Pontryagin, and it 
leads to particularly enlightening geometric formulations in many 
problems. Note that the shape of the domain of G then gives ine 
Sight into the behavior of the control program. 

In figure 1.2 a hypothetical domain of G is shown by the 
dashed line. The scalar quantity U G is maximized when G has 
a maximum projection on U at all times. Figure 1.3 shows the 
situation where more than one G vector satisfies the maximun 


condition. This generally yields a discontinuity of the control, 


hence a "corner", 


11 
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The functions B53 will be called the "forcing functions". 
It is the vector whose components are these functions which must 
be directed so as to maximize its projection on the adjoint vec- 
tor. The equation or equations which express this condition in 
the form 


Pa = h(v, Sais t) 
are obtained by ordinary calculus from (1.7). They are sometimes 


called "steering equations", 


1.7 The Hamiltonian. 
From the foregoing, one sees that the product UF has special 
significance when U is chosen according to (1.5). For convenience, 


we shall define this quantity as the Hamiltonian function, H. 
H = UF. 
The maximum principle becomes 


as) H = max , 
p 


and the general condition for an extreme of H is 
(1.9) H =0, 


which is equivalent to (1.7). 


For any one extremal, if t does not appear explicitly in the 


13 








equations of motion, the Hamiltonian is constant. This may be 


shown readily, for 


H = UF + UF 
= + + a 
UF + UFS UF P UF, 
= + oy 7: 
(U UF. )F UFP UF, . 


But U + UF, = 0 by equation (1.5), and UF, = 0 by (1.7), hence 
(le0') HU Se 


In a great many problems, the independent variable does not 
appear explicitly in the equations of motion, and the Hamiltonian 
becomes a constant of the system. In the next section it will be 
shown that in this case if the final value of the independent vari- 
able, T, is free (that is, it is not the quantity to be optimized 
and is not constrained) then H = 0 everywhere on a given path, 
provided the Euler-Lagrange equations (1.5) and (1.8) are every- 


where satisfied on the path. 


1.8 Systems which are linear in the state variables. 


Consider the functional 
f ue iy 
UF dt =fH dt =Q(T). 


If the differential equations of motion are linear in the state 


variables, then the fimctional does not involve the state variables, 
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and an absolute maximum may be obtained. If the differential 
equations are not linear, then there are no corresponding gen- 
eral proofs, and we can look only for a stationary value and a 
local extreme for the functional. 

A system of differential equations which is linear in the 
state variables with constant coefficients may be written as 

S =AS +Q 

where S is the nxl matrix of equation (1.1), Q is an nxl matrix 
whose ith component is the forcing function of the control vari-~ 
able in the ith direction, and Ais an nxn matrix of constants. 
Now, in analogy with the derivation of equation (1.6),we may form 


the integral 


i @ 
fe U(S ~ AS = Q)dt. 
O 


Integrating by parts 
a T_. 
lus} = J K+ UA )S + UQ]dt. 
O 


Taking U to be the solution to the adjoint equation U+UA = 0, we 


have the form of the Green-Lagrange formula 
T At 
Ca) fus] = f UQdt. 
oo 
Suppose that we have by some means discovered a set of 


initial values for the adjoint variables and integrated forward 


until, at some time T, 
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Hl. the curve satisfies the constraints on 
the state variables and the control vari- 
able, 

H 2. the forcing function Q is such that UPis 
maximum everywhere on the curve, and 

H 3. the terminal value of the kth component 
of the adjoint vector is unity, and all 
other components are zero at t =T. 

In this case, Faulkner (ref. 5) shows that the curve furnishes an 
absolute maximum for the kth state variable under the given initial 
values and the constraints. Akheizer (ref. 6, sect. 2.15) and 
Edelbaum (ref. 7, sects. 1.1 and 1.2) discuss conditions for 


strong and weak extremes in this and in nonlinear cases. 


1.9 The transversal condition. 
The extremals have been discussed, as the solution curves 


which satisfy the system 


(1.12) S.F=0 
U+UF =0 
S 


H =UF_ =0 
p Pp 


Consider the system already discussed, and the corresponding 


Green-Lagrange variational equation (1.6), 


T T me T 
k 
[uss | ae J (UF 6P)at = E ur] a = f $y at. 
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We require that the term on the left, evaluated at the starting 
point be zero. If s.(o) is specified, then és, (0) = Ok EES 
is not specified, then we shall choose the particular solution 
to the adjoint so that U. (0) = 0. It is necessary that the last 
term on the right be zero, else changing dt, will offer a better 


k 
trajectory. There are variations for which §P is not "small" and 


sh 

if $H dt <0, 

fe 
but for first order effects, the only condition on §S for an ex- 
tremal is that 

ué6s) =0., 

( ie 
The form of the equation does not depend on the constraints, hence 
We Will define as admissible paths those which satisfy the constraints. 
Since the constraints must be satisfied at t = T, the endpoint S(T) 
must lie on the n - k dimensional surface 

E(S,T) = 0 


Hence, for an admissible differential at the endpoint, 


= : -OS 4, °7F . /atdat = 0 
de, IZ, de, /as ds. de, /2 |. 


Ore inewa Lrixenotation, 
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(N37 db = (B,dS + ExdT), = 0 


where dé and E, are of dimension kxl, 5, is kon and d5 is nook, 
Eg may be considered to be the gradient of E in n-space. If 
E, is bordered on the right by the column Ey» this may similar- 
ly be considered to be the gredient of E in the ntl space (‘of 


the state variables and t) which we will call VE’. 


9e,/85, » « .de,/a5, de,/2t 


V = e t.] 6 


| 2¢,/25} . « .de,/25, de,,/ at 


Now if dS* is taken to be 


G 
il 


ds, 
| aT 


(1.13) may be rewritten as 
ak = (Vi*dS* ), = Q, 


In the same manner, we may have the gradient of M in ntl space, 


VM*¥ = M/asy sM/ds,, ne . d/as, nat 


and the augmented adjoint vector 


U*x = fa, Use $7 ee | 
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The vector U* must lie in the manifold spanned by the k+l vectors 


Ve,, VM , hence it may be expressed 


where Cis a lx (k+l) row vector of constants and Z* is the 
(k+1) x (ntl) matrix obtained by bordering VE* below withVM*. 
Since M is assumed to be independent, Z* has rank ktl. The 


transversal condition may be stated as the condition that the 


matrices VE*, Z* and Y* 


Y* = , a (k+2) x (ntl) matrix 
U* 
all have rank k+1. 
As a simple example, suppose we have the state variables 
W, X, y and z, with all initial conditions specified and for 
t=I w constrained, x to be maximized, y, T and z free. We may 


write 


W Be z, Jk 
ib 0 0 0 QO constraints 
0 ab 0 0 QO maximum 


Uy Uo U3 w, -H adjoint variables 


19 








a 3x5 matrix which must be of Pamke2 at t = T. Hence 


U, =u, = —H=Oatt=T. 


In the development above, free use was made of material from 

11 
a research paper of Faulkner. ‘When each e. and M involves only 
one variable of the set (S,t), the condition may be summarized as 


follows: 


il 
© 


(1.14) s.(T) free implies u, (T) 


il 
oO 


T free implies H(T) 


20 








Chapter II 
SOLUTION OF A GENERAL OPTIMUM PROGRAMMING PROBLEM USING 


DIFFERENTIAL CORRECTIONS 


2el The missing constants of the adjoint set. 


It has been shown that if relations (1.12) and (1.14) are satis- 
fied, then the resulting path generally furnishes an optimum to the 
state S(T). Finding this particular path depends upon finding a 
particular set of constants, which are the initial values of the 
unknown adjoint variables. This being true, suppose for the moment 
that the means exist for correcting a set which is "pretty good", 
until the final "perfect" set is obtained. That is, we visualize 
in the domain of U(0),the adjoint vector evaluated at t=0, the 
point U*(0) whose coordinates are these "perfect" starting values. 
If U(O) = U*, then integration of the first two of equations (1.11) 
with the control chosen according to the third of equations (1.12) 
will produce the optimum path. At some t=T on this path, the cone 
straints will be satisfied, the desired variable will be maximiz- 
ed, and simultaneously (1.12) and(.14) will be satisfied. 

In Chapter III a method will be given that will allow us to 
correct from an arbitrary point U'(0) toward U*(0) provided U'(0) 
is within some unknown region R in the vicinity of U*(0). The 
difficulty is to find some point v0) which lies within R. The so- 


called Direct Method of solution used here reduces the problem to 


all 
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the solution of a two-point boundary-value problem involving the 
extremals., Consequently, the missing boundary conditions must 
somehow be guessed. In general, the guess must be within some 
region R, if the corrective program is to converge. In the Grad- 
ient, or Steepest Ascent methods of Kelley and Bryson a nominal 
trajectory is guessed and the two-point boundary-value problem 

is solved by correcting along paths which are not extremals. This 
method too has its disadvantages. For some problems the optinum 
path may never be reached, although the correction program "con- 
verges", that is, satisfies the convergence criterion. For the 
problem treated in this thesis, a path was obtained by the direct 
methods given here which attained a final value of the variables 
to be maximized well beyond that of the "optimum" path obtained 


by the method of Steepest Ascent. 


2.2 Nominal paths. 

An advantage of the Gradient Method is the following. In 
physical optimization problems we often have some idea beforehand 
what kind of control program is likely to come fairly close to the 
optimum one. The only "guessing" involved in the Gradient Method 
is to choose such a program and construct a nominal path which 
connects the given initial constraints with the given terminal 
constraints. Provided the guess is good enough, the path is cor. 
rected by making the integrated value of 4. smaller and smaller 


along the path. 
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In a given problem, the degree of "goodness" required of this 
nominal path when using the Gradient Method corresponds in a rough 
way to the "goodness" of the first guess for U'(0) when using extre- 
mals. In the former case, however, one can be guided by physical 
reasoning. The physical meaning of the value of U'(0) is by no 
means as clear. 

In the preliminary work associated with this thesis, the authors 
have worked on a means to generate a point U'(0) within R from a 
nominal path. In the sample problem undertaken, the method works 
well and provides a value for U'(0) which is quite close to U*. 

The method is quite simple and straightforward and adds very little 
complexity to the computer program which 1s used to correct from 


Eo) to UFO). 


23 The fundamental adjoint set 
As a primary tool, we make use of the fundamental set which 


is used by Faulkner (5). The fundamental set is the nm matrix, 


@.1) 
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whose value at t=0 is the nx identity matrix and whose elements 


are chosen so that ie is the j'th component of the i'th solution 
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to the adjoint equation (1.5). Note that each row vector of the 
fundamental set is a particular solution to the adjoint (1.5). 


For instance, 


(2.2) U0) = [uy uo+-- l= for00...o]. 
0 


t= 
Now, consider a particular combination of the row vectors Uss 


(2.3) U=¢,U, te>Uet..- +e Ur 


where the c's are constants, 


This may be written as 


(2.4) U = CU, where 
55) C= fe, Co © ¢ 6 -a\: 


Note that as thus defined, U(0) = C, so that the C's are the initial 


values of the adjoint variables. At any time t, Us is the ith com- 


ponent of U: 
= + + + 
C20) B=, wy, terug +t... Fe u.,, 
or 
Gay, u, = Cu,, where 


B is the ith column of the fundamental set. 


One of the c's may be chosen arbitrarily. If some state variable 
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does not appear explicitly in the equations of motion, then the 
corresponding adjoint variable is constant. If this state vari- 
able is to be maximized then it will be convenient to make the 
corresponding c equal to plus one. If the variable is free at 
either end point then the corresponding adjoint variable is zero 
everywhere, and a different c will be chosen arbitrarily. In any 
case, we must be sure that the corresponding adjoint variable is 


noneZero. 


2.4 Program for generating the C-vector from a nominal path 

In this section a program is given for generating a starting 
set of constants U'. A nominal control program is chosen from 
physical reasoning which appears likely to be fairly close to 
optimal and which will generate a path from the given initial 
conditions to a terminal point which satisfies or nearly satis- 
fies the terminal constraints. This may be simple or difficult, 
depending on the problem. If such a program cannot be found, 
then this method for starting cannot he used, nor can be Gradient 
Method, since both are identical up to this point. In this case, 
we must go back to the "guessing game". It may be necessary to 
map" the U(0) hyperplane to find which regions give good start- 
ing points, and then to try various points within those regions. 
That is, guess the c's and calculate the corresponding trajec~ 


tories in some systematic way. 
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Suppose however that a satisfactory nominal control program 
is found. The equations of motion are now integrated numerically 
from the given initial conditions using the nominal control pro 
gram, stopping when a terminal constraint is met. At the same 
time, we integrate also the né equations of the fundamental set. 

At the terminal point S(T), each of the elements of the fundamental 
set has some value vu, 5(T). 

If the path had been the optimum path, equations (1.14) 
would have been satisfied at the point S(T). We will choose the 
vector C so that they are satisfied, using the values from the nome 
inal path. Then this vector C becomes U'(0) for starting the cor- 
rective iterations. If the program does not converge, the nominal 
path may be varied somewhat to get another trial C. In this way, 
the problem of guessing the initial values of the adjoint variables 
may be eliminated, as in the Gradient Method and still the optimum 
path is approached through extremals. 

It is to be noted that equations (1.14) are not the only 
relationships which are used to furnish the required number of 
equations in the C's. According to the particular problem, the 
Hamiltonian may be constant, or a particular adjoint variable may 
be constant. These allow use of relations calculated for t=0 


rather than at t=T as above. 
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2o5 Faulkner's method of differential corrections. 

A method was given above for obtaining a first approximation 
for the starting values of the adjoint variables from a nominal 
path. Let us now suppose that such a set is in hand which was 
obtained in this way or in some other way. 

If the equations of motion (1.1) together with the adjoint 
equations (1.5) are now integrated numerically with the control 
chosen according to the maximum principle or (1.7) a path results 
which is an extremal in that it furnishes an optimum to some end 
state. But the path obtained using this particular set of initial 
conditions does not, in general, satisfy the terminal constraints 
(1.2) and the transversal conditions (1.14). We must have some 
means for correcting these constants (which we have variously called 
the Cevector or the vector U(0) in earlier sections) so that (1.2) 
and (1.14) are satisfied at some point S(T). 

The method used here is the method of differential corrections 
developed by Faulkner (ref. 5) which makes use of differentials in 
an iterative scheme which is of the Newton-Raphson type. Making use 
of the fundamental set and the development of Chapter I, the method 
may be presented quite simply . 

Equation (1.6) provides the relation (for a problem without cor- 


ners ) 


T T 
[ss] = f UF _SPdt. 
: p 

O 0 
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For every state variable s, which is specified at t = 0, 4s,(0) = 0. 
For those which are free, we will require u;(0) = 0, so that u, 6S4 = 0 
for all cases. This procedure takes care of separated end conditions 
in a Simple and straightforward manner. This case causes a great deal 
of difficulty with some methods of solution. Note that there are 
still n unknowns at the start, since the value of s,; is unknown but 
the value of the corresponding u, is known. 


With this simplification, we have 


(2.8) [u és], = (- (UFy) SP dt, 


O 


where the Gimensions are nxn for U, nx1 for §s , nxm for Fo and mx1 
for P. From (1.7) we will choose the p's so that 

Hy = UF = 0, 
Now if S and t are fixed and the c's are varied by small amounts 5c, 


then P must vary so that 


T = 
(2.9) OcU,F, + PH, = 0, 
where the superscript T indicates the transposed matrix. But since 
U= (CU, then U, =U. Making this substitution, and noting also 


that Hop is the symmetrical mxm matrix whose ij'th element is 


=9¢ 
(Hop 43 ae) H/@p; oP, 9 


we now form the mx1 transposed of (2.9). 


Tee 2 
Zo) (UF) dc° + Hop éP = 0 
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This equation is solved for ®P, and the result is substituted 


tou 268)). 


= -1 ess 
OF = = (H) (UF) dc 


i =i r ae 
(2.11) E $s | Ri -J (iE Gi) oa dc dt, 


a form convenient for machine programming. 
Now suppose each term in the integral on the right be denoted 


as Lass Wwe have then at t =T 





U4 Sy le :° 11, dc 
ei2) ; | ; 
| nl : 5 n| Lin e t] e “nn | 5 Cn 


Since one c is arbitrary, say c,, then dc, = 0 and the equation 
furnishes n-l relations among the $§c's and the Ss's. 
In addition, the transversal conditions provide n-k equations 


whose variational form is, for the problem to be considered here, 
gu, = dcu.. 


Now for corrections we let 
a) gu, = -u,(T) - u,(T) 6r = da, 
and also, for the k constrained variables, 


és, = - 5,(2) - 5,(t) $1 + sQ(7) 
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where S067) is the given terminal constraint for the ith variable. 
Since there are k of these variables constrained, there remain 
(n.k-1) elements bs, in equation (2.12). These may be eliminated, 
yielding k equations. The n-k equations (2.13) are adjoined to 
these, for a total of n equations involving the un!mowns OT, oC, 
bc, . . . 68, in terms of the integral elements I, 5. the 
elements of U(T), the deriviatives at t=T and the differences (2.13). 


The equations are solved for the unknowns, and these are applied as 


corrections for the start of the next iteration. 
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Chapter III 


DEVELOPMENT OF THE PROBLEM 


3.1 General description of the problem and assumptions. 

The problem described in this thesis is that of deter-= 
mining the optimal angle of attack program p(t), which will 
maximize the distance covered over the earth's surface by a 
hypersonic glider or lifting vehicle which has been injected 
MmeavQ sOme initial re-entry point. The starting potmt vor ec 
trajectory is specified by the re-entry injection parameters, 


mee initial conditions, which are: 


v(0) = 25,920 feet/second - initial velocity 
h(0) = 30€ ,000 feet - initial altitude 
= x(0} = O nautical miles - initial distance 
z(O0) = 0.18 degrees - initial flight path angle 


These parameters were taken from reference (8) so that a 
comparison could be made between the Gradient Method used by 
Carlier authors and the method of differential corrections 
which was used in this thesis. Also, the value of wingload= 
ing, the value of acceleration due to gravity, the value for 
a standard earth radius, and the value of a standard nautical 
mile were taken from the same reference. 

Wingloading = mg, = 27.3 1b/ft* , where: 


A 


m mass of the vehicle 


wing plan form area 


A 
go = bec ft/sec” - acceleration due to gravity at 


the earth's surface 
The model atmosphere used was based on the tables of reference 


(9). The atmospheric density was assumed to have the form 
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p=ne where the parameter b was calculated to fit the tab- 
ular values at h = 0 feet and h = 200,000 feet. It was fur- 
ther assumed that the atmosphere was spherically symmetric and 
fixed with respect to the earth for simplifying reasons. 
fee Hovations of motion. 

the coordinate system uséd 15 G@eeiepecmeae woe ede ae 
equations of motion are: 


x VECOsmr., 


Bs 

R+h 
: n 

7 ) 


Vo Sa Zen 


<je 
i 


ie 
-~7—- g cos Zz, 


ta 


v g 
— Ne —_ Z, 
vb mV ata Bal COS Z@ v cos 


Where: R = 3440 nautical miles, tne radius of the earth, 


1 
D=slpv acy), 
eee 
L=3l(pv AC), 

R42 
B= Eola 


The coefficients of lift and drag are from reference (8) and : 


are: 
Cy = Cro Sine COs ge sin p | 
ae = 6 lsinpl+ G 
a 9} ae ne ; DO 

Where: Cr, = 1.82 ; Cy, = Wao: Cha = 0,042 ; and p is the 


angle of attack of the vehicle. 


3.3 Adjoint equations. 
‘he system of differential equations which is adjoint 


to the variational equations of the equations of motion (3.2) is: 
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ew eae 1D _ é& 
up, = bG@Esne cos z|u,, + l= ah on 222 zu, 
rn E COS.Z 2) 08 @ as 2282) y 
(R+h)- mv oh ch V Zz 
« [= 126682 de 
me) 6, [so |v, -|sin z|U, - E elu, 
L Civ]. . COsez een cos |Z 
7 Ex: ~ Ov mv oR +h Vv Ju, 
‘aac Rv sin ay Sin 2) u al " 7 
7 = | a ” | cos Z| pai le cos Z| uy, 
Vern Zoe 2 sin 2 
4 | R+ hn V Ju, 


The fundamental set of adjoint equations U or equivalently 
(u, ;) was chosen such that the fundamental set at time t = 0 


Memeene identity matrix of rank four. That ic: 


u(y =f 2 2 of 
ee A ow 

(3.4) Lu 2 3 a 
UE We Uy 





. De ee 7 : 
and further: (uy 4) 16 =f, (u, ,)dt , where the Ws are defined 
by the adjoint equations (3.3). We have the further relation 
Aor 2 rit 
that: U= [u, up u, ui] = [e, c, C. c,, ][U]” which is a more 
convenient form for computations using fortran programming. 
Tn general, the solution of the adjoint eauations is 
only determined to within a multiplicative constant , hence 
we can choose one of the c's. Tn this problem, we chose C,=1. 
Tf we consider the equations of motion in vector or 


matrix form, we have the following relationships: 





X : 0 
=" h = — = @) = 
S — i = = i 
ve as BD ok where F. = & 
Z =m p cp 
Lp 
mv 
The Hamiltonian is given by H = UF, 


3.4 Maximum principle. 

For any extremal, we must satisfy the condition that 
UF = 0 for any point on the extremal. This condition 
uniquely determines the angle of attack program for that 


particular extremal since: 


ad CL eae) 
al cry pallet, 4 ae 

Usk, = apuz VSpUy = 0 

Substituting for = and 7 and solving for p we obtain: 
2 
i) Pe eaGeren sat v | 1° 8| “10, age 
Ze u_ Wy a 9iC Z, V 
LO “Zz DL 





W 
fOr = a < JO 


aoe) ena conditions. 

Since we wish to maximize the total distance traveled 
over the earth's surface, we must have the following conditions 
meeevme cna time t = 

ey) maximum 
(4720 } h(T) 


v(T), 2z(T), and T free. 


0 


I 


3.6 Variational eauetions. 
On any curve without corners, we have the following 
relation between the adjoint and the variations: 


m 
he 


J J J J - As: _ 
307) [us &x + uy, Oh + usev + ud6zJn= 4 "op dtiwejieaive. 3 


For an extremal we must also satisfy the Euler equations or 
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the maximum principle of Pontryagin so that we have 

ew od = 
fur dt = maximum and U:F = O along the extremal. If we 
0 


consider x,h,v,z, and t fixed anc change the constants then 


—_ — => 
CU OU 0U i a 
Ee OC, + %e.,°°3 + @ toy| 7 + Usk Cp = 702 Oring ther 


° 3 oc Dp pp 
eu _ pi ay | | 
de, = , 1 =2,3,4, We can solve for Gp and subsureuLre in 
equation (3.7) which will give us four equations with 


which we can correct the initial values of the c's for a 


terminal variation of the altitude h: 


Lb Tae) 1 
j j | j j ' . = he aa 
(3,8) judex + uy oh + vam He upoz aa a, (U Ba) i atée, 
UF op 
T (GF) (WE) 
Let 1,,-/ ——P.—2atte 
0 U-F 
pp 
= 1 a 5 h 
since uy = ilies ee O for all t, and noting that 


Th; = Bay 
in terms of the oon 


ve can solve the following set of equations for 6h 


2 2 Z 
~ 3 caer Pee 
Li 
Phy ae OZ Tou Toy Tyudp [ecu 


From the transyersal conditions, we have that: 
uy (T) = QO 
(3.10) u,(T) = 0 
(Xu thu, ) ap = 0 
From the last equation of (3.10) we have: uy, (T) = (-cot Z) on 
=a : 
Since the Hamiltonian H = U°F and H(T) = 0 from equations 
(2710). then the Hamiltoniianmiemidenvieally zeno for alla. 
This fact will be used to determine one of the constants (Co) 
and effectively reduce the problem of guessing the initial 


Gs) 60 start the calculations. 
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The transversal conditions (3.10) give us three more 
independent variational equations associated with the end 


Ceomed tions:: 


wr 
a 
] 


2 Bice | 
[uvéc, a wybe, 5 fen Jin 


c 


Ey In 


ee 
be 
j- 
a a) 
S 
a 
{| 
NG < & 


[uséc, + uzbe, + U 
é[ xu,+ oni = h(T)Lurec, + upbe. a: usbey Je 
3.7 Differential corrections 
The Wethod of “differential corrections 2-1 e, 
Maulkner in reference (5) was used for the solution of the 
Sovemization problem. In general, we will euess the™inieia) 
values associated with the adjoint enuations and the terminal 
time T. For this problem, we must choose Co, C3: Cy, ‘nd T. 
We then correct these varameters by calculating am extremal 
fame che initial guesspandc vtvhenssolwve for corrections oe 
parameters in terms of the end values associated with the 
extremal. For all end conditions of the form s(T) =S, we will 
generally have some error in s(T). We set 6s(T) = S = s(T) = 
s(T) 6T. For this problem, we set 6hn(T) = n(T)-h(T)éT and in 


a similar manner, we obtain from equation (3.11): 


2 3 4 ? cn \ 
Lule. & u,éc. a Uy bey Ip By be CE ee 


t 
Ss 
| 
I 


, Fema 5 by ° 
Pere) «6©-u_(T) = Luréc, + uzeo, - uC) Jin as u(T)ET 


hy 


r a a fod 3¢ 
cot Z Jip = h(T) upto, sp up fe. ey 
+( xu, + hu, + hay | pET 
These equations together with the solution to (3.9) for €h in 


terms of the &c's will give a set of equations from which we 
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solve for the ¢c's and €T to make corrections in the starting 


values for the subsequent trajectory. 
3.8 Computation of approximate startine values by use of 
nominal trajectories. 

In general jy ror a nth” ordertsev or Oli ferential 
equations, one must choose (n-1) of the constants of the sol- 
ution to the adjoint equations to start the voroblem. This is 
Bpenerally a difficult problem in itself, since the constants 
are not related to the ohysical aspects of the problem in any 
discernable manner. One of the primary aspects of this thesis 
WaeeecvO find a method which could be used to find a logical 
Choice for the initial valves of these constants. Fortunately 
Pie Hamiltonian was identically zero in this problem. Tnis 
fact, together with the transversal conditions (3.10) enabled 
mer cO sOlve for approximate constants Dy first compuvi ise 
nominal trajectory; here we choose a@ nominal angle of attack 
program p{t), which seemgslikely in a physical sense to 
give a miximun distance trajectory. The first such nominal angle 
of attack program tried was to use the maximum lift over drag 
ratio where p=20.5 degrees, a const-unt. This did not work. 

We then used the following angle of attack program for which 
we are indebted to Professor Faulkner: 


ieee Sur [2 ~ h(t) | 


750,000 
The value of 57” for p is the value at which we have 
maxium lift. It was suvposed that the argle of attack was 
never negative and from the data of constant angle of attack 
programs, the value of 750,000feet was chosen to insure that 


p was never negative. 
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solve for the ¢c's and €T to make corrections in the starting 


values for the subsequent trajectory. 
3.8 Computation of approximate startine velues by use of 
nominal trajectories. 

In general, for a nth order set of aifterential 
equations, one must choose (n-1) of the constants of the sol- 
ution to the adjoint equations to start the voroblem. This is 
generally a difficult problem in itself, since the constants 
are not related to the vohysical aspects of the problem in any 
discernable manner. One of the primary aspects of this thesis 
was to find a method which could be used to find a logical 
choice for the initial valves of these constants. Fortunately 
Me Hamiltonian was identically zero in this proolem:. aes 
fact, together with the transversal conditions (3.10) enabled 
Hee tO SOlve for approximate constants by Tirsv Compuvraees 
nominal trajectory; here we choose a nominal angle of attack 
program p(t), which seemslikely in a physical sense to 
Bae a maximum distance trajectory. The first such Nomina ane. 
of attack program tried was to use the maximum lift over drag 
ratio where p=20.5 degrees, a constant. This did not work. 

We then used the following angle of attack program for which 
we are indebted to Professor Faulkner: 


einai anit ” Be 80 | 


The value of Tc for p is the value at which we have 
maxium lift. It was suvposed that the argle of attack was 
never negative and from the data of constant angle of attack 


programs, the value of 750,000feet wus chosen to insure that 


p was never negative. 
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With this choice for the angle of attack program, we 
obtained a good set of starting constants by integrating the 
fundamental set U along with the ejuations of motion until h = O. 


At this time, we solve the set of equations: 


2 3 ly ji 
Uy Uy UD, Co “Uy, “cot Z 
2 3 Ly 7 L 
13) uy uy, wy e3| = [-uy 
2 3. O44 a 
U u U Cc -U 
Z Z Z Im Ly Zz ~ 


ome CD , C3; and cy. These esuations are the transversal 
equations @.10). This problem of guessing the c's is further 
reduced due to the Hamiltonian being identically zero, since 
mom a given choice of C3 and cy, the angle of attack is umigely 
determined from equation (3.5 and then we may solve for co from 
the Hamiltonian at time t = 0. If we have the case where n(O)= 
O, then we have a singularity and connot solve for Co, DUG 
jas case, the choice of co is arbitrary. Bitectavely ae. 
problem is then reduced by an order of one and we need only 
guess two constants. We felt that if we had an extremal, we 
woulda determine its constants in this manner and ence we cour 
get estimates for the constants from a ee a "good" 
aporoxination to an extremal. 

3.9 Computational scheme. 

This problem was vrogrammed on a CDC 1604 computer 
using a Runge-Kutta integration method. The comoutational 
scheme is as follows: 

i) We need a nominal trajectory to start. To get it, we 
first intergrate the equations of motion@2) and the fundamental 


set (3.4) using a programmed angle of attack and ending the 
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trajectory when h = 0. At this time, we solve for the 
initial constants cp, cz, and cy, from equation (3.13) We then 
vse these constants to start the first optimal trajectory. 

ii) Then we compute the first optimal trajectory where 
the angle of attack program is determined by the maximum 
principle. We integrate the equations of motion, the fundament= 
al set, and the nine Tay in equation (3.9) a total of twenty= 
five equations. We terminated the integration when one of 
the following conditions is met: h = 0, uu, = OF wu = oor 
lzlz us Since it was felt that the end of the steric ite... 
of the trajectory had been reached. We felt that in this 
problem NOLOPTIMUM. VOULd Ceccure: oF uL< O since the resulting 
fen drag would dissioate energy needlessly. The supposition 


that this is a maximum distance trajectory dictates that l2i< & 


© 


rl 


Since we do not have t appearing explicitly in the equat= 
ions and since the terminal conditions determine T, we do not 
use €T in our calculations. ‘We only need the 6c's the correct 
the c's to start the next trajectory. 

411) We then repeat the computations as in the above 
paragraph, using the new c's for the next trajectory. ke 
moo form convergence at the end Of “eacen trajec or,. 
(3.14) te (h(t)}* + [ku, + hula + fuj(7)]* + fun)? €¢, 

we say that we have converged to the solution. 
3.10 Results and conclusions. 

The optimum trajectory that was calculated is depicted 
on Fig.33 and the corresoonding angle of attack program is 
denicted on Fig. 3.4. On these two figures, a comparison has 


been made between the results of this thesis and results from 


reference (8), which were obtained by the Gradient (or Steev- 
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est Ascent) Method. The comparison can rot be considered as 
an exect one since the integration method, the time step 
or integration interval, and the apnvroximation to the 
atmospheric models may not be identical. However, there 
appears to be a significant increase in the maximum distance 
obtained using the method of differential corrections. 

The convergence criteron G@ 14) we initially chose turned 
out to be a poor choice. We could not determine how close we 
were to the optimal solution. We then proceeded to "map" 


the plane of c3 and cy by using a systematic choice of c3 and 


Cy amc computing the cormmesponding extvremals. Ttemresm a a 
the mapping procedure are given on Fig.3.5, Fig.3. ,and Fig.3.7- 
Fig.3.5 shows the contours of distances obtained by extremals 
as a function of C3 and cy). Fig.36 gives detailed contours 
Gm@eund the optimal solution. Fig.3.7 is a cross section of ine 
Prirape” which occurs in the C4 and c, Dlane. It is a con- 
jecture of the autnors that some difficulties are encountered 
by a gradient method solution when there exists such a ridge 


in aeproolem, since the slope is nearly zero. 


The major problem encountered was the apoarent instability 
of the vroblem near the end of an extremal. We attempted to 
integrate backwards from the terminal conditons and found that 
this was virtually impossible. The problem seems to be unstable; 
we could not even integrate backwards with constant p. Back- 
ward integration was very good if we started backward from a 


meme thet was not too close to the end of thewura jccv. 


Another main problem encountered was that of determining 


LO 





the end time on 4 trajectory. We first thought that all we 
need consider is the time when h = 0 . However, we found 
tna other conditions effectively terminated the useful part 
of the extremal earlier. For example, most extremals led to 
a point where lzl> = 

Theoretically, on the maximum distance trajectory, all 
terminal conditions would be met simultaneously. In the actual 
computations, this did not occur. The most vrobable reason 
that we have this apparent discrepancy is that the round off 
errors in the calculations and the inherent error in the 
integration method are large enough to prevent us from satis= 
fees 2ll conditions at the end time, 

In the final analysis, we found that a program which 
meanuced the magnitude of the corrections would lead to con 
wereence. The method used was to calculate 4 set of Cc ™s 
from a nominal trajectory and then calculate the first tra- 
jectory and the ¢c's corresnronding to that trajectory. We 
store x, the c's, and éc's of this trejectory and then cal-— 
culate the next trajectory. If the distance of this trajectory 
is greater than the distance of the preceding one, we proceed 
with the iterative process. If the distence is less, then 
we reduce the corrections by one half and calculate the tra- 
Weoceory azain. if we still have Vess ais7aic eee 1 Ur Uaer 
reduce the corrections until we obtain a greater distance. 

The reason we developed this method was because the magnitudes 
of the éc's were of the order of magnitudes of the c's [ft 
was noted that the integrals relating Gh to the Ges are 


improper intezrals. This, in addition to the unstable char- 
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acter of the end of the trajectory, tended to make the 
magnitude of the tc's large; however, the program did 
correct in the right direction. These integrals seem to be 
divergent for the extremals which satisfy the transversal con- 
ditions. A formulation of the problem suggested by Professor 
Paulkner which would avoid these divergent integrals by using 
a set of u's which were determined at the end time was pro- 
eeemmed., However, this program was not chéekedgwout or run 
saumee the original integrals did not actuslly diverac lime i 
culations made on the computer. The method sufszestea was to 
wise a fundamental set U such that U (T) = Cag, After com> 
puting an extremal using U , we solve for the matrix Of Comes 
stants Lay 5) which satifies the relation: u(T) = 4, ,JU(T). 
We then compute the same extremal using ifs vu. + ct. 


i 1 pug 
U, => a.,U’, i= 1,2. We then have the relation 
1 j=] ij 


where 


Typ" = \2 
Eh(T) = if (Uy ae dt &c and we solve for éc to correct the 
-__--— 


O59 + "Pop 
starting vaiuves. Further, this integral is not divergent. 

The first integration method used in celculations was 
rectangular integration anc it was found thet trajectories 
thus calculated were very rough. A comvarison of the rectang- 
ular integretion method and the Runge-Kutta method was mede 


ana we chose the latter method due to its accurscy and 


smoothness. 
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